# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse", "estimatr", "reshape2"))

library(tidyverse)
library(estimatr)
library(reshape2)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

long_mturk <- 
  mturk_opeds %>%
  filter(responded_w3 & responded_w2) %>%
  select(subject_id, Z, contains("scale")) %>% 
  select(-contains("mj")) %>%
  melt(id.vars = c("subject_id", "Z")) %>%
  mutate(dv = sapply(strsplit(x = as.character(variable), split = "_"), function(x)x[2]),
         wave = sapply(strsplit(x = as.character(variable), split = "_"), function(x)x[4]),
         study = "mturk")

long_elite <- 
  elite_opeds %>%
  filter(responded_w2) %>%
  select(subject_id, Z, contains("scale")) %>% 
  select(-contains("mj")) %>%
  melt(id.vars = c("subject_id", "Z")) %>%
  mutate(dv = sapply(strsplit(x = as.character(variable), split = "_"), function(x)x[2]),
         wave = sapply(strsplit(x = as.character(variable), split = "_"), function(x)x[4]),
         study = "elite")

long_df <- bind_rows(long_mturk, long_elite)

gg_df <- 
  long_df %>%
  group_by(Z, dv, wave, study) %>%
  do(tidy(lm_robust(value ~ 1, data = .))) %>%
  ungroup %>%
  mutate(Z = factor(Z, levels = c("control", "amtrak","veterans", "climate","wallstreet","paul"),
                    labels = c("Control", "Treatment: Amtrak", "Treatment: Veterans", 
                               "Treatment: Climate", "Treatment: Wall Street", "Treatment: Flat Tax")),
         dv = factor(dv, levels = c("amtrak","vets", "wall","flat", "climate"),
                     labels = c("Outcomes: Amtrak", "Outcomes: Veterans", 
                                "Outcomes: Wall Street", "Outcomes: Flat Tax",
                                "Outcomes: Climate")),
         wave = as.numeric(as.character(factor(wave, levels = c("w1", "w2", "w3"), labels = c("0","10", "30")))),
         study = factor(study, levels = c("mturk", "elite"), labels = c("MTurk Study", "Elite Study")))


g <- 
  ggplot(gg_df, aes(x = wave, y = estimate, group = Z, color = Z, fill = Z, shape = Z)) +
  geom_point(size = 4) +
  geom_line(alpha = .1) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha =.2, show.legend = FALSE) + 
  theme_bw() +
  facet_grid(dv ~ study) +
  ylab("Composite Attitude Scale") +
  xlab("Days Since Treatment") +
  guides(linetype = "none", fill = "none", alpha = "none", ncol = 2) +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_blank())
g

# Show that W1 effects don’t vary by reporting status ---------------------

fit_1 <- lm_robust(dv_amtrak_main_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_2 <- lm_robust(dv_climate_main_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_3 <- lm_robust(dv_flat_main_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_4 <- lm_robust(dv_vets_main_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_5 <- lm_robust(dv_wall_main_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))

fit_6 <- lm_robust(dv_amtrak_scale_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_7 <- lm_robust(dv_climate_scale_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_8 <- lm_robust(dv_flat_scale_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_9 <- lm_robust(dv_vets_scale_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))
fit_10 <- lm_robust(dv_wall_scale_w1 ~ Z*responded_w2, data = filter(mturk_opeds, responded_w1))

